\documentclass[12pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath, amssymb}
\usepackage{geometry}
\usepackage{lmodern} % 使用更现代的字体
\geometry{a4paper, margin=1in} % 设置页面边距

% 中文支持
\usepackage{ctex}

% 定义一些命令
\newcommand{\controlinput}[1]{\ensuremath{\mathtt{#1}}}
\newcommand{\statevar}[1]{\ensuremath{\mathit{#1}}}
\newcommand{\outputvar}[1]{\ensuremath{\mathbf{#1}}}

\begin{document}

\title{基于《质子交换膜燃料电池系统建模与仿真》的功率控制方程组推导}
\author{物理学教授（讲解原理）与自动化教授（讲解控制）}
\date{}
\maketitle

\section*{引言：系统的“自动化”视角}

\textbf{自动化教授：} 各位同学，今天我们的课题是如何为一套质子交换膜燃料电池（PEMFC）系统建立一个用于控制的数学模型。我们的目标非常明确：我们希望通过调节三个我们能直接操纵的“控制输入（Inputs）”，来精确控制系统的最终“输出（Output）”——也就是净功率。

\begin{itemize}
    \item \textbf{控制输入 (U)}:
    \begin{itemize}
        \item 氢气入口摩尔流量 ($\controlinput{q_{H_2}^{in}}$) —— 我们的主要“油门”。
        \item 空气入口摩尔流量 ($\controlinput{q_{air}}$) —— 我们的“节气门”，必须与油门匹配。
        \item 冷却水摩尔流量 ($\controlinput{q_W}$) —— 我们的“散热器”控制器。
    \end{itemize}
    \item \textbf{控制输出 (Y)}:
    \begin{itemize}
        \item 系统净输出功率 ($\outputvar{P_{system}}$) —— 这是客户真正得到的功率。
    \end{itemize}
\end{itemize}

我们的任务就是建立连接 $U$ 和 $Y$ 之间的动态方程组链条。这个链条中会有一个核心的“状态变量（State）”——电堆温度 $\statevar{T}$，它会产生系统的主要动态（延迟）响应。

---

\section*{链条一：辅助功耗（寄生负载）模型}

\textbf{物理学教授：} 首先，我们必须清醒地认识到，这个系统不是一个“永动机”。为了让它“呼吸”，我们必须输入空气，而这需要空压机，它会消耗我们自己发的电。这在热力学上称为“寄生负载”。

\textbf{自动化教授：} 对。我们的控制输入 $\controlinput{q_{air}}$（空气流量）不是凭空而来的，它需要消耗功率。我们需要先对这个功耗建模。

\subsubsection*{1. 压缩机绝热功率 $P_{ad}$}

\textbf{物理学教授：} 根据热力学定律，将气体从入口压力 $p_s$ 压缩到出口压力 $p_d$，在绝热过程中所需的功率由下式给出：
\begin{equation*}
P_{ad} = \frac{\kappa}{\kappa-1} p_s q_v \left[ \left( \frac{p_d}{p_s} \right)^{\frac{\kappa-1}{\kappa}} - 1 \right] \quad \text{(基于公式 1)}
\end{equation*}
\begin{itemize}
    \item \textbf{现实含义}：这个公式计算了理想情况下压缩空气所做的功。空气流量越大（$q_v$），或者压缩比（$p_d/p_s$）越高，所需的功率就越大。
    \item \textbf{参数/变量}：
    \begin{itemize}
            \item $p_s$：入口压力（通常是大气压）。
            \item $p_d$：出口压力（即阴极入口压力 $p_c$）。
            \item $q_v$：空气的体积流量 ($m^3/s$)。我们的控制输入是摩尔流量 $\controlinput{q_{air}}$，它们可通过气体状态方程 $q_v = \controlinput{q_{air}} \cdot (RT/p_s)$ 转换。
            \item $\kappa$：空气的绝热指数（约为1.4）。
        \end{itemize}
\end{itemize}

\subsubsection*{2. 压缩机实际功耗 $P_{comp}$}

\textbf{自动化教授：} 理想总是美好的，但现实中没有100\%的效率。我们必须考虑压缩机的“绝热效率” $\eta_{ad}$。
\begin{equation*}
P_{comp} = P_{ad} / \eta_{ad} \quad \text{(公式 2)}
\end{equation*}
\begin{itemize}
    \item \textbf{现实含义}：这是我们为了维持 $\controlinput{q_{air}}$ 的空气流量，实际需要从电堆“偷”走的电功率。
    \item \textbf{参数}：$\eta_{ad}$ 是绝热效率，一个0到1之间的值（本文中取0.75-0.85）。
\end{itemize}

---

\section*{链条二：电化学主功率模型}

\textbf{自动化教授：} 这是我们的核心。我们将展示如何从“油门” $\controlinput{q_{H_2}^{in}}$ 一步步推导出电堆的“毛功率” $P_{stack}$。

\subsubsection*{1. 从氢气流量到电流 $I$}

\textbf{物理学教授：} 根据法拉第电解定律，电流的产生本质上是电子的定向移动，而电子的来源是氢气的电化学反应（$2H_2 \rightarrow 4H^+ + 4e^-$）。因此，电流和消耗的氢气量是严格成正比的。
\begin{equation*}
I = \frac{2 \cdot u \cdot F}{N} \cdot \controlinput{q_{H_2}^{in}} \quad \text{(基于公式 3)}
\end{equation*}
\begin{itemize}
    \item \textbf{现实含义}：这是模型的关键入口点！它将我们的控制输入 $\controlinput{q_{H_2}^{in}}$ 直接转换为了系统的内部电流 $\statevar{I}$。
    \item \textbf{参数}：
        \begin{itemize}
            \item $u$：氢气利用率（一个效率系数，表示有多少比例的氢气被成功反应）。
            \item $F$：法拉第常数（$96484600 \, C/kmol$），代表每摩尔电子的电荷量。
            \item $N$：电堆中串联的单电池数量（本文为332）。
        \end{itemize}
\end{itemize}

\subsubsection*{2. 内部气体分压 $p_{H_2}, p_{O_2}$}

\textbf{物理学教授：} 有了电流 $\statevar{I}$ 和我们输入的空气流量 $\controlinput{q_{air}}$，我们就可以计算出在稳态下，电极表面（催化剂层）的反应物“浓度”——也就是气体分压。

\textbf{自动化教授：} 请注意，这篇论文的模型在这里做了简化。它使用的是代数方程，而不是微分方程。这意味着模型假设气体分压能够瞬时响应电流和温度的变化，而没有模拟气体的传输延迟。

\begin{enumerate}
    \item \textbf{水蒸气饱和压力 $p_{H_2O}$} (只与温度有关):
    \begin{equation*}
    \log(p_{H_2O}) = 2.95\times10^{-2}(T-273.15) - 9.18\times10^{-5}(T-273.15)^2 + 1.44\times10^{-7}(T-273.15)^3 - 2.18 \quad \text{(公式 4)}
    \end{equation*}

    \item \textbf{氢气分压 $p_{H_2}$} (阳极):
    \begin{equation*}
    p_{H_2} = 0.5 p_{H_2O} \left[ \frac{1}{\exp\left(\frac{1.635(I/A)}{T^{1.334}}\right) \frac{p_{H_2O}}{p_a}} - 1 \right] \quad \text{(公式 5)}
    \end{equation*}

    \item \textbf{氧气分压 $p_{O_2}$} (阴极，方法一):
    \begin{equation*}
    p_{O_2} = p_{H_2O} \left[ \frac{1}{\exp\left(\frac{4.192(I/A)}{T^{1.334}}\right) \frac{p_{H_2O}}{p_c}} - 1 \right] \quad \text{(公式 6)}
    \end{equation*}
\end{enumerate}
\begin{itemize}
    \item \textbf{现实含义}：这些经验公式描述了一个核心现象：电流 $\statevar{I}$ 越大，反应越剧烈，电极表面的反应物分压（$p_{H_2}, p_{O_2}$）就越低。同时，温度 $\statevar{T}$ 越高，水蒸气压力 $p_{H_2O}$ 越大，也会影响反应物分压。
    \item \textbf{参数/变量}：$T$ 是电堆温度（来自热力学模型），$A$ 是电极面积， $p_a, p_c$ 是阳极和阴极的入口总压（由 $\controlinput{q_{air}}$ 等因素决定）。
\end{itemize}

\subsubsection*{3. 氧气浓度 $C_{O_2}$}

\textbf{物理学教授：} 活化损失（下一步）与“浓度”更相关。根据亨利定律，溶解在电解质中的氧气浓度 $C_{O_2}$ 与其气体分压 $p_{O_2}$ 成正比。
\begin{equation*}
C_{O_2} = p_{O_2} / ( (5.08 \times 10^6) \exp(-498 / T) ) \quad \text{(公式 11)}
\end{equation*}

\subsubsection*{4. 单电池电压 $U$ (电压损失模型)}

\textbf{物理学教授：} 电池的实际输出电压 $U$ 等于其理论最高电压 $E_{Nernst}$ 减去所有内部损失。
\begin{equation*}
U = E_{Nernst} - U_{act} - U_{ohmic} - U_{con} \quad \text{(公式 14)}
\end{equation*}
\textbf{自动化教授：} 这是关键的一步，它告诉我们电压 $U$ 是如何被电流 $I$ 和温度 $T$ 共同决定的。

\begin{itemize}
    \item \textbf{理论电压 $E_{Nernst}$} (Nernst Voltage):
    \begin{equation*}
    E_{Nernst} = 1.229 - 0.85\times10^{-3}(T-298.15) + 4.3085\times10^{-5}T[\ln(p_{H_2} \cdot p_{O_2}^{0.5})] \quad \text{(公式 9)}
    \end{equation*}
    \item \textbf{活化损失 $U_{act}$} (Activation Loss):
    \begin{equation*}
    U_{act} = -[-0.9 + 0.0024 T + 7.3\times10^{-5} T \ln(C_{O_2}) - 0.00023 T \ln(I)] \quad \text{(公式 10)}
    \end{equation*}
    \item \textbf{欧姆损失 $U_{ohmic}$} (Ohmic Loss):
    \begin{equation*}
    U_{ohmic} = (0.01605 - 3.5\times10^{-5} T + 8\times10^{-5} I) \cdot l \quad \text{(基于公式 12)}
    \end{equation*}
    \item \textbf{浓差损失 $U_{con}$} (Concentration Loss):
    \begin{equation*}
    U_{con} = -B \ln\left[1 - \frac{J}{J_{max}}\right] \quad \text{(公式 13)}
    \end{equation*}
\end{itemize}
\begin{itemize}
    \item \textbf{现实含义}：
        \begin{itemize}
            \item $E_{Nernst}$ 是理想电压，由反应物分压（$p_{H_2}, p_{O_2}$）和温度 $T$ 决定。
            \item $U_{act}$ 是启动化学反应所需的能量损失，受 $T$, $C_{O_2}$ 和 $I$ 影响。
            \item $U_{ohmic}$ 是电流 $I$ 流过电堆内阻（受 $T$ 影响）时的电压降。
            \item $U_{con}$ 是在高电流 $J$（$J=I/A$）下，反应物来不及扩散到催化剂表面造成的电压骤降。
        \end{itemize}
\end{itemize}

\subsubsection*{5. 电堆总功率 $P_{stack}$}

\textbf{自动化教授：} 有了单电池电压 $U$ 和总电流 $I$，计算总功率就水到渠成了。
\begin{align*}
U_{stack} &= N \cdot U & \text{(公式 15)} \\
\outputvar{P_{stack}} &= U_{stack} \cdot I & \text{(公式 16)}
\end{align*}
\begin{itemize}
    \item \textbf{现实含义}：这是电堆产生的毛功率。
\end{itemize}

---

\section*{链条三：热力学动态（核心反馈回路）}

\textbf{物理学教授：} 我们不能忘记能量守恒！电堆产生的总能量，一部分以电能 $P_{stack}$ 输出，剩下的大部分都转化为了废热 $\Delta Q_{stack}$。这些热量必须被有效管理，否则电堆会烧毁。
\begin{equation*}
\Delta Q_{stack} = N \cdot I \cdot (1.481 - U) \quad \text{(公式 17)}
\end{equation*}
\begin{itemize}
    \item \textbf{现实含义}：电流 $I$ 越大，或者电压 $U$（效率）越低，产生的废热就越多。
\end{itemize}

\textbf{自动化教授：} 这里就是我们系统的动态核心，也是我们的第三个控制输入 $\controlinput{q_W}$ 发挥作用的地方。我们将建立一个热平衡方程，它在Simulink中通常被实现为一个积分器，从而引入延迟响应。

\subsubsection*{动态热平衡方程}

\textbf{物理学教授：} 系统温度 $\statevar{T}$ 的变化率 ($\frac{dT}{dt}$) 取决于“净热流”。
\begin{equation*}
C_{th} \frac{d\statevar{T}}{dt} = \dot{Q}_{in} - \dot{Q}_{out}
\end{equation*}
\textbf{自动化教授：} 基于论文中的公式(17)到(20)，我们可以构建这个动态方程：
\begin{center}
    \textbf{（标注：以下为基于原文公式推导出的动态微分方程，原文(20)式是其稳态形式）}
\end{center}
\begin{equation*}
C_{th} \frac{d\statevar{T}}{dt} = \underbrace{(\Delta Q_{stack} + Q_{gas\_in})}_\text{热量来源} - \underbrace{(\Delta Q_{rad} + \Delta Q_W + Q_{gas\_out})}_\text{热量散失}
\end{equation*}
\begin{itemize}
    \item \textbf{热量来源}：
        \begin{itemize}
            \item $\Delta Q_{stack}$：电堆产生的废热（公式 17）。
            \item $Q_{gas\_in}$：随入口气体带入的热量（通常较小）。
        \end{itemize}
    \item \textbf{热量散失}：
        \begin{itemize}
            \item $\Delta Q_{rad}$：通过辐射散失到环境的热量（公式 18）：
            \begin{equation*}
            \Delta Q_{rad} = \delta \cdot \sigma_b \cdot a \cdot (T^4 - T_0^4)
            \end{equation*}
            \item $\Delta Q_W$：被冷却水带走的热量（由 $\controlinput{q_W}$ 控制）（公式 19）：
            \begin{equation*}
            \Delta Q_W = \controlinput{q_W} C_{p,w} (T - T_0)
            \end{equation*}
            \item $Q_{gas\_out}$：随出口气体带走的热量。
        \end{itemize}
    \item \textbf{参数}：$C_{th}$ 是系统的总热容（一个表征热惯性的关键参数，在动态模型中必须存在）。$\delta, \sigma_b, a, C_{p,w}$ 都是热力学相关参数（发射率、玻尔兹曼常数、面积、水比热容）。
\end{itemize}
\textbf{自动化教授：} 这个微分方程就是我们系统的“状态方程”。它告诉我们，当前的温度 $\statevar{T}$ 是过去所有热量输入和输出累积（积分）的结果。这就产生了热延迟！当电流 $I$ 突然增加时，$\Delta Q_{stack}$ 瞬时增加，但 $\statevar{T}$ 只会缓慢上升。

---

\section*{最终输出：系统净功率 $P_{system}$}

\textbf{自动化教授：} 最后，客户得到的净功率 $\outputvar{P_{system}}$ 是电堆产生的毛功率 $\outputvar{P_{stack}}$ 减去我们为了维持系统运行而消耗的寄生功率 $\outputvar{P_{comp}}$。
\begin{equation*}
\outputvar{P_{system}} = \outputvar{P_{stack}} - \outputvar{P_{comp}} \quad \text{(基于图4)}
\end{equation*}
\begin{itemize}
    \item \textbf{现实含义}：这就是我们的最终控制目标。
    \item \textbf{控制哲学}：我们的控制器必须同时调整 $\controlinput{q_{H_2}^{in}}$ 和 $\controlinput{q_{air}}$。如果 $\controlinput{q_{air}}$ 太低，会导致“氧气饥饿”，$P_{stack}$ 会崩溃。如果 $\controlinput{q_{air}}$ 太高，$P_{comp}$ 会剧增，导致 $P_{system}$ 反而下降（如论文图4所示）。同时，控制器还必须调整 $\controlinput{q_W}$ 来稳定 $\statevar{T}$，因为 $\statevar{T}$ 会反馈影响链条中的几乎每一个方程。
\end{itemize}

\end{document}
